Adapted from: Floudas1999; Section 3.5 and Lasserre2009; Table 5.1
Consider the polynomial optimization problem from Floudas1999; Section 3.5
A = [
0 0 1
0 -1 0
-2 1 -1
]
bz = [3, 0, -4] - [0, -1, -6]
y = [1.5, -0.5, -5]
using DynamicPolynomials
@polyvar x[1:3]
p = -2x[1] + x[2] - x[3]
using SumOfSquares
e = A * x - y
f = e'e - bz'bz / 4
K = @set sum(x) <= 4 && 3x[2] + x[3] <= 6 && f >= 0 && 0 <= x[1] && x[1] <= 2 && 0 <= x[2] && 0 <= x[3] && x[3] <= 3
Basic semialgebraic Set defined by no equality 8 inequalities 4.0 - x[3] - x[2] - x[1] ≥ 0 6.0 - x[3] - 3.0*x[2] ≥ 0 24.0 - 13.0*x[3] + 9.0*x[2] - 20.0*x[1] + 2.0*x[3]^2 - 2.0*x[2]*x[3] + 2.0*x[2]^2 + 4.0*x[1]*x[3] - 4.0*x[1]*x[2] + 4.0*x[1]^2 ≥ 0 x[1] ≥ 0 2.0 - x[1] ≥ 0 x[2] ≥ 0 x[3] ≥ 0 3.0 - x[3] ≥ 0
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.Optimizer
Clarabel.MOIwrapper.Optimizer
A Sum-of-Squares certificate that p≥α over the domain S, ensures that α is a lower bound to the polynomial optimization problem.
The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
end
solve (generic function with 1 method)
The first level of the hierarchy gives a lower bound of `-7``
model2 = solve(2)
nothing # hide
-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 9
constraints = 12
nnz(P) = 0
nnz(A) = 24
cones (total) = 2
: Zero = 1, numel = 4
: Nonnegative = 1, numel = 8
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 2.2646e+00 2.2646e+00 0.00e+00 4.87e-01 2.46e-01 1.00e+00 2.84e+00 ------
1 4.8146e+00 4.8773e+00 1.30e-02 8.53e-02 3.55e-02 2.33e-01 6.39e-01 8.00e-01
2 5.8894e+00 5.8986e+00 1.56e-03 1.07e-02 4.45e-03 3.33e-02 9.41e-02 8.57e-01
3 5.9970e+00 5.9972e+00 3.18e-05 2.02e-04 8.47e-05 6.56e-04 1.82e-03 9.81e-01
4 6.0000e+00 6.0000e+00 3.18e-07 2.02e-06 8.47e-07 6.56e-06 1.82e-05 9.90e-01
5 6.0000e+00 6.0000e+00 3.18e-09 2.02e-08 8.47e-09 6.56e-08 1.82e-07 9.90e-01
6 6.0000e+00 6.0000e+00 3.18e-11 2.02e-10 8.47e-11 6.56e-10 1.82e-09 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 24.9ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -6.00000e+00
Dual objective value : -6.00000e+00
* Work counters
Solve time (sec) : 2.49300e-02
Barrier iterations : 6
The second level improves the lower bound
model4 = solve(4)
nothing # hide
-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 82
constraints = 101
nnz(P) = 0
nnz(A) = 242
cones (total) = 10
: Zero = 1, numel = 20
: Nonnegative = 1, numel = 1
: PSDTriangle = 8, numel = (10,10,10,10,...,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 9.7795e-01 9.7795e-01 1.11e-16 6.56e-01 4.75e-01 1.00e+00 2.05e+00 ------
1 1.2254e+00 1.2350e+00 7.85e-03 1.39e-01 9.77e-02 1.75e-01 4.77e-01 7.96e-01
2 2.6674e+00 2.7439e+00 2.87e-02 6.47e-02 4.85e-02 1.71e-01 2.05e-01 7.36e-01
3 4.3661e+00 4.4172e+00 1.17e-02 1.41e-02 9.22e-03 7.83e-02 4.58e-02 8.30e-01
4 5.1102e+00 5.1418e+00 6.18e-03 7.04e-03 4.30e-03 4.73e-02 2.38e-02 5.82e-01
5 5.5386e+00 5.5472e+00 1.54e-03 2.08e-03 1.20e-03 1.33e-02 7.23e-03 7.49e-01
6 5.6457e+00 5.6476e+00 3.40e-04 4.70e-04 2.64e-04 3.00e-03 1.66e-03 8.07e-01
7 5.6906e+00 5.6907e+00 1.34e-05 1.82e-05 1.02e-05 1.19e-04 6.47e-05 9.72e-01
8 5.6922e+00 5.6922e+00 8.64e-07 1.17e-06 6.55e-07 7.63e-06 4.16e-06 9.40e-01
9 5.6923e+00 5.6923e+00 1.21e-08 1.62e-08 9.08e-09 1.06e-07 5.77e-08 9.87e-01
10 5.6923e+00 5.6923e+00 2.66e-10 3.50e-10 1.97e-10 2.33e-09 1.25e-09 9.78e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 301ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -5.69231e+00
Dual objective value : -5.69231e+00
* Work counters
Solve time (sec) : 3.00959e-01
Barrier iterations : 10
The third level improves it even further
model6 = solve(6)
nothing # hide
-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 451
constraints = 506
nnz(P) = 0
nnz(A) = 1376
cones (total) = 10
: Zero = 1, numel = 56
: PSDTriangle = 9, numel = (55,55,10,55,...,55)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 6.2703e-01 6.2703e-01 0.00e+00 7.75e-01 5.19e-01 1.00e+00 1.48e+00 ------
1 7.4953e-01 7.6135e-01 1.18e-02 1.77e-01 1.32e-01 2.39e-01 4.26e-01 8.06e-01
2 1.2333e+00 1.3121e+00 6.39e-02 4.83e-02 5.13e-02 1.73e-01 1.70e-01 8.68e-01
3 1.9671e+00 2.0062e+00 1.99e-02 1.35e-02 1.14e-02 6.32e-02 3.87e-02 8.18e-01
4 2.7047e+00 2.7536e+00 1.81e-02 9.27e-03 5.40e-03 6.99e-02 1.94e-02 7.43e-01
5 3.0495e+00 3.0614e+00 3.90e-03 4.17e-03 2.24e-03 2.13e-02 8.43e-03 9.22e-01
6 3.4706e+00 3.4774e+00 1.96e-03 1.08e-03 4.82e-04 9.63e-03 1.95e-03 8.04e-01
7 3.7081e+00 3.7129e+00 1.30e-03 5.58e-04 2.22e-04 6.51e-03 9.26e-04 6.34e-01
8 3.9158e+00 3.9173e+00 3.71e-04 1.47e-04 5.33e-05 1.93e-03 2.38e-04 8.24e-01
9 3.9469e+00 3.9488e+00 4.92e-04 1.00e-04 2.67e-05 2.34e-03 1.30e-04 6.86e-01
10 4.0073e+00 4.0086e+00 3.31e-04 4.62e-05 8.58e-06 1.52e-03 4.86e-05 7.76e-01
11 4.0267e+00 4.0275e+00 1.76e-04 2.25e-05 4.45e-06 8.15e-04 2.55e-05 7.35e-01
12 4.0582e+00 4.0584e+00 3.91e-05 4.13e-06 7.46e-07 1.79e-04 4.47e-06 9.70e-01
13 4.0660e+00 4.0660e+00 1.14e-05 1.12e-06 2.01e-07 5.20e-05 1.22e-06 8.03e-01
14 4.0671e+00 4.0671e+00 6.05e-06 5.62e-07 1.01e-07 2.75e-05 6.16e-07 6.69e-01
15 4.0683e+00 4.0683e+00 7.91e-07 7.24e-08 1.31e-08 3.60e-06 7.98e-08 8.78e-01
16 4.0684e+00 4.0684e+00 3.27e-07 2.87e-08 5.20e-09 1.48e-06 3.17e-08 7.67e-01
17 4.0685e+00 4.0685e+00 6.90e-09 6.03e-10 1.09e-10 3.13e-08 6.67e-10 9.80e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 57.0ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -4.06848e+00
Dual objective value : -4.06848e+00
* Work counters
Solve time (sec) : 5.69970e-02
Barrier iterations : 17
The fourth level finds the optimal objective value as lower bound.
model8 = solve(8)
nothing # hide
-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 1736
constraints = 1855
nnz(P) = 0
nnz(A) = 5436
cones (total) = 10
: Zero = 1, numel = 120
: PSDTriangle = 9, numel = (210,210,55,210,...,210)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 5.8462e-01 5.8462e-01 1.44e-15 8.33e-01 6.41e-01 1.00e+00 1.32e+00 ------
1 6.6260e-01 6.6529e-01 2.69e-03 1.88e-01 1.72e-01 2.55e-01 4.10e-01 7.95e-01
2 8.8639e-01 9.8255e-01 9.62e-02 7.03e-02 8.83e-02 2.42e-01 2.31e-01 7.79e-01
3 1.3550e+00 1.3811e+00 1.92e-02 1.49e-02 1.66e-02 5.46e-02 4.98e-02 8.09e-01
4 1.9495e+00 1.9837e+00 1.75e-02 8.55e-03 7.91e-03 5.63e-02 2.40e-02 7.05e-01
5 2.2256e+00 2.2451e+00 8.77e-03 6.02e-03 4.95e-03 3.46e-02 1.55e-02 5.32e-01
6 2.7604e+00 2.7709e+00 3.82e-03 1.85e-03 1.19e-03 1.56e-02 4.09e-03 7.86e-01
7 2.9706e+00 2.9778e+00 2.42e-03 1.20e-03 6.94e-04 1.07e-02 2.53e-03 4.62e-01
8 3.2375e+00 3.2402e+00 8.29e-04 3.95e-04 1.91e-04 3.90e-03 7.88e-04 8.75e-01
9 3.2912e+00 3.2953e+00 1.26e-03 3.05e-04 1.22e-04 5.47e-03 5.55e-04 5.46e-01
10 3.3766e+00 3.3796e+00 8.92e-04 2.36e-04 9.20e-05 4.08e-03 4.36e-04 4.69e-01
11 3.6119e+00 3.6139e+00 5.40e-04 8.55e-05 2.73e-05 2.43e-03 1.46e-04 7.14e-01
12 3.7063e+00 3.7081e+00 4.86e-04 5.81e-05 1.67e-05 2.18e-03 9.52e-05 6.34e-01
13 3.9300e+00 3.9304e+00 1.02e-04 1.29e-05 3.31e-06 4.89e-04 2.08e-05 8.27e-01
14 3.9879e+00 3.9880e+00 1.64e-05 2.15e-06 5.26e-07 8.02e-05 3.48e-06 8.64e-01
15 3.9979e+00 3.9979e+00 3.11e-06 3.51e-07 8.41e-08 1.48e-05 5.72e-07 9.11e-01
16 3.9996e+00 3.9996e+00 7.01e-07 7.14e-08 1.71e-08 3.29e-06 1.18e-07 8.93e-01
17 3.9999e+00 3.9999e+00 1.22e-07 1.23e-08 2.94e-09 5.71e-07 2.03e-08 8.38e-01
18 4.0000e+00 4.0000e+00 2.04e-08 1.59e-08 4.74e-10 9.50e-08 3.25e-09 9.18e-01
19 4.0000e+00 4.0000e+00 3.69e-09 2.84e-09 8.55e-11 1.72e-08 5.87e-10 8.30e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 807ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -4.00000e+00
Dual objective value : -4.00000e+00
* Work counters
Solve time (sec) : 8.06567e-01
Barrier iterations : 19
This notebook was generated using Literate.jl.